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SUMMARY 


The classical theory of flight dynamics fcr airplane longitudinal 
stability and control analysis was extended to the case of a hypervelocity 
reentry vehicle. This Induced the elements Inherent In supersonic and 
hypersonic flight such as the Influence of the Mach nimwer on aeroc^namlc 
characteristics, the effect of reaction control system and aerodynamic 
controls on the trim condition through a wide range of the speed. 

Phugold motion and angle-of-attack oscillation for three typical cases 
of cruising flight, ballistic entry and glide entry v«re Investigated. In 
each case, closed form solutions for the variations In altitude, flight path 
angle, speed and angle-of-attack were obtained. These solutions displayed 
explicitly the influence of different design parameters and trajectory 
variables on the stability of the irotlon. 

Aerodynamic and Inertia characteristics of a Single-Stage-To-Orbit 
vehicle (SSTO) were used In the numerical application but the basic analytical 
theory developed can be used to assess entry dynamics of any lifting vehicle. 
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I. INTRODUCTION 


A new Earth orbital transportation system could be required beginning In 
the last decade of this century provided that such system offers clear and 
significant cost/performance advantages over the existing system. 

A proposal for such system Is a S1ngle-$tage-To-0fb1t vehicle (SSTO) 
capable of performing space missions similar to that of the current space 
shuttle. After a direct ascent Into orbit and completion of Its space mission, 
the vehicle must sustain a reentry Into the Earth's atmosphere at hypersonic 
speeds with lifting capability for a crossrange of 1100 n.ml., and perform 
a horizontal landing at speed not exceeding 170 knots. 

Studies associated with the entry dynamics of this vehicle are now being 
conducted. Several excellent analyses have demonstrated the applicability 
of the control configured design approach to Advanced Earth Orbital Trans- 
portation Systems. It has been shown that, fOr the control -con figured 
vehicle (CCV), the major design criteria are the hypersonic trim, the location 
of the center-of-gravlty, the control surfaces heating and the mass of the 
system. Applicability and potential performance gains of CCV design concepts 
as applied to the development of a SSTO were globally analyzed In reference 1. 
One of the most significant applications of CCV concept Is the reduction of 
the required vehicle static margins which can be used to reduce stabilizer 
and fin size, reduce trim drag, reduce wing size because of more favorable 
balancing tall loads. This may lead to significant weight savings. This 
problem has been analyzed by Freeman and Wilhite (ref. 2) using modern design 
analysis techniques (refs. 3-4). Starting with a baseline venicle with 
positive level of longitudinal stability, the vehicle was configured with 
the conventional stability requirements relaxed. The control -con figured 
vehicle has an estimated gross lift-off mass 10 percent less than that of 
the baseline vehicle design. The largest single component weight reduction 
Is a 33 percent reduction In the wing weight due to the 43 percent reduction 
In exposed wing area. Because of this, for the CCV, the ratio of body- 
pi anform to wing-pl anform area Increases and the aerodynamic center Is much 
more dependent on Hach nimber. From the aerodynamic characteristics deduced. 

It has been shown that the control -configured vehicle can be trimmed hyper- 
sonlcally over an acceptable angle-of-attack range, from 10® to 60®, by a 
combination of a body flap and elevon deflection. 

On the other hand, for the SSTO vehicle, an aft center-of-gravlty loca- 
tion resulting from large engine mass at the rear of the vehicle has created 
new stability problems during entry. The effect of center-of-gravlty loca- 
tions on hypersonic trim has been analyzed in reference 5. The results of 
this excellent study predict that control configuring of the reentry vehicle 
Is definitely a requirement, and a combination of aerodynamic controls and 
reaction control system (RCS) would be adequate. While the location of the 
center-of-gravlty can be easily evaluated, and by averaging, the RCS charac- 
teristics can be accurately modeled, the aerodynamic characteristics, and In 
particular, the aerodynamic center, can only be projected from wind tunnel 
tests before definite results can be obtained from flight tests data. This 
results in uncertainties to be Incorporated In any analysis of dynamic 
stability and control of advanced transportation systems, lliese are analyzed 


for the space shuttle orbiter in references 6 and 7 and can be extended to 
the SSTO design. 

In this report* we shall extend the classical theory of flight dynamics 
for airplane stability and control analysis to the case of a SSTO vehicle. 

This will Include the elements Inherent In supersonic and hypersonic flight 
such as the Influence of the Mach number on aerodynamic characteristics* 
the effects of RCS and aerodynamic controls on the trim condition through 
a wide range of the speed. Another Important aspect of stability and control 
In hypersonic flight such as In the case of a SSTO vehicle which Is not 
present In the case of airplane dynamics Is the strong Influence of the 
trajectory variables. While at low speeds one can always assess the stability 
dynamics using small perturbations theory about a steady reference flight 
path* reentry dynamics Is an unsteady phenomenon In which the effects of the 
deceleration* the atmospheric density gradient and the curvature of the 
Earth are all Important. 

The organization of this report Is as follows. After this Introductory 
chapter* the classical Euler's equations of motion for airplane stability 
analysis are derived In Chapter II. For entry dynamics* these equations 
must be modified to include the effect of the variation of the vector 
acceleration of the gravity. As a first step* we Introduce the modifying 
term to obtain the equations for longitudinal dynamics to be used for the 
analysis In this report. The effectiveness of reaction control system and 
aerodynamic controls on the trim condition during ^ntry Is analyzed In 
Chapter III. Entry dynamics of the SSTO Is the subject of Investigation In 
the following two chapters. In Chapter IV the phugold* or trajectory 
oscillations for the three typical cases of cruising flight* ballistic entry 
and glide entry are Investigated. In the analysis the Important effects of 
atmospheric gradient* curvature of the flight path and deceleration are 
displayed explicitly. Chapter V Investigates the angle-of-attack oscilla- 
tions during entry for these cases. Resonance phenomena are analyzed and 
criteria fOr dynamic stability are established. 

This study benefits from the wealth of aerodynamic data for the control - 
configured SSTO vehicle generated In references 2 and 8 and because of that 
It provides further understanding of the dynamic behavior of this advanced 
transportation system although the basic analytical theory can be used to 
assess entry dynamics of any lifting vehicle. 
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II. EULER'S EQUATIONS OF MOTION 


The best description of the dynamic forces and nwmen^is connected with 
a rotating rigid body Is given by Euler's dynamical equations based upon an 
axis system that remains r1g1d1y fixed at the center-of> gravity of the moving 
body. (N» the other hand, the notion of the center-of-gravlty Is generally 
described by a differential system written with respect to Inertial system. 

In this chapter, we shall discuss the relationship between the two 
systems of equations. 

Equations of Motion Over a Flat Earth 

For a non rotating fiat Earth, the variables Involved are: 

a) Position of the center-of-gravlty with respect to a Cartesian 
ground system GxgygZg considered as Inertial system 



Notice that H Is the altitude and subscript g h^s been omitted for convenience, 

b) Components of the velocity vector ^ along the body axes 




u 

V 

w 


c) Euler's angles 
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d) Components of angular velocity H of the vehicle along the tK)dy axes 


P 


H 


q 

r 


3 


With respect to the body axis system, the velocity vector V is defined 
in direction by the angle-of-attack a and the sideslip angle B with the 
relations (Fig. 1) 

u = V cosB cosa 

V » V sinB 

w = V cosB sina 

V * u cosa COSB + V sinB + w sina cosB (1) 

= / + v^ + w^ 


tana = w/u 
tanB = v//u^ + w^ 

where V is the^magnitude of the velocity vector and u, v and w are the 
components of V along the body axes. 

Using the body axes as reference system, the two vector equations 


♦ i5x?) = t 
^ M = I 


( 2 ) 


for the forces and moments have projections 
m[u + qw - rv] = X - mg sine 

m[v + ru - pw] = Y + mg cose sin(|i (3) 

m[w + pv - qu] = Z + mg cose coS(|) 

and 

IxP - Ixz(?‘ + Pq) + dz - = Mx 

lyq + dx ^ ^xz(P^ “ 

Izf - Ixz(P - rq) + (ly - Ix)pq = 

In Eqs.(3), X, Y and Z are the components along body axes of the combined 
aerodynamic and propulsive forces. In Eq.(4), it is assumed that, by 

symmetry the products of inertia Ixy = lyz = 0- The three moments of force 

My, My and are about the roll, pitch and yaw axes respectively. 

We have, in addition, the kinematic relations 


4 


t i 



Fig. 1. Body axis system. 
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( 5 ) 


* (q sin 4 t + r cos<t>)/cose 

6 =« q cosij) - r sinij) 

I = p + (q s1n<|) + r cos<j))tane 

The body system CXYZ Is obtained from the ground system 6xyz by three 
rotations: about z to Cx^y^z, then e about yi and then 4 > about X, 

(Fig. 2). Hence, we have the transformation 


X ^ 


’ cosi|) -si nip 0 ’ 


cose 0 sine 


10 O' 


■ X ■ 

y 

= 

s1nip cosip 0 


0 1 0 


0 coS(p -si nip 


Y 

z 

N J 


0 0 1 

k J 


-sine 0 cose 


0 si nip cosip 


Z 

L / 


Carrying out the matrix multiplication 


X 1 


cosecosip (sinesinipcosip-cosipsinip) 

(slnecosipcosip+slnipsinip) 

r 

y 

= 

cosesimp (si nesinipslmp+cosipcosip) 

(si necosipsimp-sinipcosip) 


. z . 


-sine cosesinip 

cosecosip ^ 

V 


Putting (X,Y,Z) = (u,v,w) and (x,y,z) = (ug.vg.wq), we deduce the projections 
of the velocity vector on the ground Inertial axis system 

Xg = I’g = u cosecosij/ + v(s1nes1n(j(Cosi|; - cosij)Slnii») 

+ w(s1n0cos<j)cosij^ + s1n(|)sintfi) 

yg = Vg = u cos0s1ni(( + v(s1n0s1n4»s1nij; + cos(|>cost;;) (6) 

+ w(s1n0cos<|)s1ni|; - s1n(|)Cos4;) 

Zg = -ft = Wg = -u s1n0 + V cos0s1n(|) + w cos0cos(j) 


The equations (3), (4), (5) and (6) completely describe the motion of 
the vehicle In the case of airplane dynamics. 

The aerodynamic and propulsive forces and moments are functions of the 
control Input, the speed and altitude. On the right-hand-side, we have 9 
variables 


H, u, V, w, 0, (j), p, q, r 

It suffices to integrate the equations concerning these variables. Once they 
are known as functions of the time, the additional equations In x, y and 
can be Integrated. 

Longitudinal Equations of Motion Over a Spherical Earth 

For airplane dynamics, the equations above are adequate. For any vehicle 
configuration with known mass and engine properties, It suffices to tabulate 
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the aerodynamic characteristics for a simulation study to assess the dynamic 
response under any flight condition. This approach Is Illustrated In 
reference 9 where stall and post-stall characteristics of a high-performance 
fighter airplane, the F-16, with relaxed longitudinal static stability were 
analyzed. 

For long-range hypervelocity vehicle, as In the case of the present 
SSTO vehicle, a flat Earth formulation Is clearly Inadequate because of the 
fact that the gravity acceleration changes Its magnitude, and more Importantly, 
Its direction along the flight path. 

We assume that the Earth Is sphetical but non-rotating and consider the 
fight In the plane of a great circle (Fig. 3). We first start with the 
equations for flight over a flat Earth by considering the relevant variables 


w, e and q. 


m[u + qwj = X - mg sine 

(7) 

m[w - qu] = Z + mg cose 

(8) 

lyq = My 

(9) 

9 = q 

(10) 

ft = u sine - w cose 

(11) 


The original inertial ground system GxgZg is now moved by translation to the 
geocentric inertial system Oxz with the origin at the center of the Earth. 
Let R be the position vector of the center-of-gravity with respect to the 
system Oxz. For airplane motion, the equations have been derived on the 
basis that the vector g is always directed along the ground axis Szg. In 
reality, this vector has rotated by an angle x equal to the variation in 
the longitude. Let y be the angle between the local horizontal and the 
velocity vector. This angle is termed the flight path angle. We have the 
kinematic relations 

dR/dt = V siny (12) 

dx/dt = {V/R)cosy (13) 

Since the position vector rotates with angular speed x, the angular rate e 
for the rotation of the body axes, referred to the fixed axes is 

9 “ q + (V/R)cosy (14) 

This equation replaces Eq.(lO) and reduces to it in the case of flat Earth, 

R ->• «. Since ft ■ ft, by comparing Eqs.(ll) and (12), we have 

V siny ■ u sine - w cose (15) 

From Eqs.(l) for the case of zero sideslip angle, we have trivially 

u » V cosa 

(16) 

w ■ V sina 
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To this, we add the geometric relation 

0 • Y + a (17) 

By taking the derivative of Eq.(15), using Eq.(14) 

V siny + yV cosy • u sine - w cose + (u cose + w sine )0 
• (u sine “ w cose) + (u cose w sine)(q + cosy) 

Now combining Eqs.(7) and (8) 

m[(u sine - w cose) + q(u cose + w sine)] » X sine - Z cose - mg 
Using Eq.(18) we write this equation as 

m[V^ siny + yV cosy - (u cose + w sine)^ cosy] = X sine - Z cose - mg 
But 

u cose + w sine » V cose cosa + V sine sina 
* V cos (e-a) » V cosy 

Hence 

m[V siny + yV cosy - j cos^y] = X sine - Z cose - mg (19) 

On the other hand, by eliminating q between the Eqs.(7) and (8) 
m[uu + ww] » Xu + Zw - mgu sine + mgw cose 
Therefore, since = u^ + w^ and V^' » uu + ww, we have 
mVV “ XV cosa + ZV sina - mgV(sin0 cosa - sina cose) 
or 

mV » X cosa + Z sina - mg siny (20) 

By observing that 

Fj ■ X cosa + Z sina (21) 

is the projection of the combined aerodynamic and propulsive force on the 
tangent to the flight patn, we have the force equation along the tangent to 
the flight path 

mV * Fj - mg siny (22) 

Furthermore, using Eq.(20) in Eq.(19) and solving for y 
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Hence 


yl 

IlhyV COSY ■ COS*Y + X Sln0 - Z COS8 


mg 


- X coso slny - Z si no siny mg sin^y 


mVy ■ X si no - Z coso - m(g - ^ )cosy 
Observing that 

* X sino - Z coso ^23) 

is the projection of the combined aerodynamic and propulsive force along the 
normal to the flight path, we have 

mVy * ■ *"(9 “ ^ )cosy (24) 

In summary, for longitudinal motion of a hypervelocity vehicle, the 
variables are conveniently selected as 


R» V, Y, g, 0 
with equations 


f - » slnv 


dt 


m 


- g siny 


V 




yz 

(g - P )cosy 


(25) 



dt 


m 


q 



cosy 


The auxiliary variables are o, u, and w with the relations 


0 » Y + o 
u • V cosot 
w » V si not 
tana » w/u 


( 26 ) 
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III. CONTROLS EFFECTIVENESS 


In the pitching moment equation In system (25), the total pitching moment 
Is represented by 

My • H + Mt - If (I^ - Ij)s1n2e (27) 

where M Is the aerodynamic moment, Is the pitching moment due to reaction 
control system and the last term expresses the gravity torque. A complete 
discussion of gravity torques has been given by Roberson (ref. 10) and 
Beletskll (ref. 11). Because R Is large, the effect of this term Is small 
whenever aerodynamic and thrusting moments are effective. But the term should 
be retained in the case of orbital flight when M -»• 0. 

Effectiveness of Reaction Control System 

The thrusting moment Is an arbitrary function of the time. If f(t) Is 
the thrust magnitude applied at a distance from the center-of- gravity, for 
the duration t of the firing It can be represented by an average constant 
value (Fig. 4) 


a f (t) 


M. ^ 28 ) 

L To 


^ t 

Fig. 4. Average thrusters force. 



This average value accounts for the effects of thrust built up with time and 
thrust loss due to back pressure Increases with decreasing altitude. 

The total pitching moment coefficient Is defined as 

M. M. 3g(I - I )s1n2e 

^"’tot * ^ pSV^c ' ^m H pSV^c ■ R^SV^c 

where C„, Is the aerodynamic moment coefficient and the reference length c 
used is the rean aerodynamic chord. To trim at Cmtot * 

thrusters (RCS) Is effective at high altitude, while at high dynamic pressure, 
aerodynamic moment Is the main Ingredient for control. Let f be the average 
constant thrusters force. The thrust moment coefficient Is 


12 


(ycHf/is psv^) 


(30) 


t 



M, 


t 


h pSV'c 


Introducing the load factor 
, pSV^C, qSC, 

" • w “-iw ir* 


where q is the dynamic pressure, we write Eq.(30) as 


(31) 


» (yc)(f/W)(CL/n) (32) 

t 

This formula makes explicit the effectiveness of reaction control system. 

On the right-hand-side of the equation, the first ratio 1s obviously a design 
parameter while the second term, which Is simply the thrust-to-welght ratio, 
concerns the sizing of the engines. The last term depends on the trajectory 
to be flown. As a rule. In an optimum fashion, jet firings can be used as 
an effective means of pitch control whenever the load factor, or proportionally 
the dynamic pressure, remains small. For a reentry trajectory, RCS jets can 
be used to supplement pitch control during the early part of the entry at 
low dynamic pressure, but under severe trim condition their use may also be 
mandatory at lower altitude with significant penalty In fuel consumption. 

Also for a skip trajectory, at high load factor, the use of RCS requires 
high thrusting force If aerodynamic pitch control Is constrained by 
aerodynamic heating consideration. 

Effect of the Aerodynamic Center 

For the aerodynamic moment, consider a wing-body combination (Fig. 5) 
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The force system acting In symmetric flight can be represented as a 
lift L and a drag D acting at a reference point, the mean aerodynamic center, 
together with a pitching moment Mac> Then, the aerodynamic moment about the 
vehicle center-of- gravity Is 

M + (L coso + D s1na)(h - h,^)c + (L sina - D coso)zc (33) 

aC aC 


where h and h.^t expressed In fraction of the mean aerodynamic chord c. 
Indicate the location of the center-of-gravlty and the aerodynamic center, 
respectively. If the last term which Is small Is neglected, the aerodynamic 
pitching moment coefficient Is 


C 


m 



+ (h - h 


ac 



(34) 


where 


C|^ ■ Cj^ cosa + Cp si no 


(35) 


Is the normal aerodynamic force coefficient. 

The Inviscid theory of thin wing at small o predicts that Is 
Invariant with o, but In general It varies with both the angle-of-attack 
and the Mach number. Likewise, the location of the aerodynamic center 
which Is here expressed as hac In fraction of the nwan aerodynamic chord Is 
a function of the angle-of-attack and the Mach number. It is particularly 
sensitive to the variation of the Mach number In the case of the CCV 
because of the large ratio of the body-planform to wing-planform area. 

For small angle-of-attack, and considering that Coa « Cl, we have 
the classical formula for use In the case of airplane dynamics 


+ (h - h^^)C, 
m m ac L 

aC 


(36) 


Hence, on the basis that Is Independent of the angle-of-attack, and In 
the case where hac Is also insensitive to variation In a we have the formula 



(h - h 


ac 



(37) 


which provides an excellent way of finding hac from test results, that Is 
from measurements of the stability derivatives Cnu and CIq. Since CLo > 0 
and for static stability, Cm^ < 0, the quantity (n - haclt 1r fraction of 
the mean aerodynamic chord c. Is meaningfully called the stability static 
margin. 

Since In the case of hypersonic reentry. In the early part of the 
trajectory, the angle-of-attack Is large, we retain the more rigorous 
expression (34) and measure the stability static margin according to 


K • h - h » 3C_/3C» 
ac m N 


( 38 ) 
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since Cri Is a positive and Increasing function of a and at the trim condition 
^ * Cnu < 0 provides the condition for the vehicle to be statically 

stable, a positive K denotes a level of static Instability. An aft shift of 
the center- 0 f-gravlty, or during an entry a shift forward of the aerodynamic 
center will Increase this level. 


Equation (37) Is only appropriate for the computation of the value of 
K at low angle-of-attack and low flach number when K is Insensitive to the 
variations of these variables. In general, from Eq.(34', we have the equation 


K • 



7T ■ 


(39) 


which shows that, at a given Mach number. K Is a function of the angle-of- 
attack through the functions Cm(a), and CN(a). These functions can be 
modeled based on wind tunnel test data but the small and nearly constant 
value Cffljc can only be estimated with certain degree of arbitrariness. As 
a first approximation, we can neglect this term and take 

K « K • C /C„ (40) 

m N 

Typical plots of K and K are shown In Fig. 6, 


K » h-h„c 



Fig. 6. Maxlmwn shift of the aerodynamic center at 
hypersonic speeds. 
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The difference K - K ■ Cmac/C(^ becomes small at high ang1e-of*attack 
with Urge values of Cn. It Is seen that fbr both K and ? there exists 
angU'Of-attacic providing the maximum shift forward of the aerodynamic center 
yielding the maximum level of static margin. This leads to the fbllowing 
procedure for evaluating Cmgc* basis of the Invariance of Cniac A''*! 

by taking the successive partial derivatives of the equation 


Cm • ♦ KC„ 

m m. ^ n 

ac 


with respect to o wt have 

^m * '^a^N * 
a a 


and 


■ K C» + 2K + KC», 
m aa N a N N 
aa a Gux 


(41) 


(42) 


(43) 


At the point a where K Is maximized* • 0* and we have 


K 



(44) 


and hence, combining this equation with Eq.(43) 


C C 
ra m 
a _ gg 

r” z — 

''N 

g gg 




gg 


(«) 


At the point a. where K Is maximized we have 




(46) 


The partial derivatives can be evaluated If Cnj(g) and CN(g) are accurately 
modeled as continuous functions of g. E'lujtion (46), up?**! solving, provides 
g« but not the value g where K Is max1m-i7,e<1. Since the difference Cnj^^/C^ 
between the two curves Ic and K Is small, we can assume that the 
behaviors of the tw curves, respectively near the two stationary values 
are the same, that Is they have the same slope and the same curvature. 

Since Ka(a«) ■ Kjj(g) ■ 0, we have the condition for the same curvature 


K 

gg 


(g) 


gg 



Vn 


gg 


'n 


g*g© 


(47) 


Upon substituting this approximation In Eq.(45) we have the equation for 
evaluating g where K Is maximized 
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( 48 ) 






S 



Once this value of a has been computed, we obtain from Eqs.(41) and (44) the 
constant value 



'm- 


■m 

a 


(49) 


For a typical reentry profile, a constant and high angle-of-attack is 
maintained through the hypersonic range. Since at hyperst^ic speeds, the 
aerodynamic characteristics are nearlv indepen<fent of the Mach number the 
value of K is nearly constant (ref. S) and can be kept at a low level of 
instability by selecting a high enough angle>of-attack (Fig. 6). Below 
Mach 5.0, the aerodynamic characteristics are strongly influenced by the 
Mach number. This, coupled with the rapid decrease in the angle>of-attack 
for a transition frm entry to glide configuration may create unsteady 
phenomena resulting in large static instability. lT)e faction K is now 
function of the two variables o and M. Since it is not convenient to 
analyze a three-dimensional plot we consider the intersections of the surface 
K * K(a, M) by the planes a ■ constant. 

Aerodynamic data for the example CCV are taken from references 2 and 8. 
The reference dimensions are given in Fig. 7 with the center-of- gravity 
located at 0.69 of the body length. Aerodynamic forces and mcxnent are 
controlled by a combination of a body flap and an elevon system. With 
■ 0, and - 0 and for different values of the angle-of-attack, the 
variations of C|_ and Cj>| as functions of the Mach manber are shown in Fig. 8. 
At low angle-of-attack. Cm can be approximated by Ci but at high angle-of- 
attack especially in the vicinity of sonic speeds, due to a sharp rise in 
the drag, the two curves can be quite distinct. Hence, we continue to 
evaluate K through the use of Eq.(39) but now with an estimate value for 
based on the ass(am)tion that at each Mach number, it is independent of 
the angle-of-attack. Hence, we can take Cm.. • C,n 

evaluation has to be used with care because^the point corresponds to nearly 
zero angle-of-attack where the aerodynamic center is not well defined. The 
value for Cm,, obtained is compared with the value of es titrated by using 
Eq.(49) at low angle-of-attack (a ** 2.5“) and this provioes excellent agree- 
ment. The result is plotted in Fig. 9 in dashed line as function of the 
Mach number. In the same figure, we have plotted the variation of C„, for 
different values of tha angle-of-attack. From the data in Figs. 8 and 9, 
the location of the aerodynamic center, now expressed in percentage of the 
body length is evaluated through the formula 


X 


ac 



i 


(C, 


"ac ■ 




(50) 
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Mach number 


Fig. 8. Lift and normal force coefficients as 
functions of the Mach number for 6- * 0. 6. * 0 




Mach number 


Fig. 9. Aerodynamic pitching moment coefficient 
as function of the Mach number. 





with a xcg/t * 0.69, the location of the aerodynamic center as function of the 

Mach number is shown In Fig. 10 for low angle-of-attack and In Fig. 11 for 

high angle-of-attack. These figures make explicit the Influence of the 
Mach number on the longitudinal stabiMty characteristic. 

During entry, the vehicle Is slightly unstable at hypersonic speeds. 

It can even be made slightly stable with the selection of high angle-of- 
attack for entry. Below Mach 5.0, the aerodynamic center shifts forward 
and the vehicle becomes more unstable. Then at lower Mach number, the aero- 
dynamic center begins to shift aft and the neutral point, defined as ® Xcg, 
is reached at a Mach nuirtber slightly above 2.0. For all the angles-of- 
attack considered, the maximum level of stability which corresponds to the 
maximum aft-location of the aerodynamic center occurs in the transonic 

region. This is due to a sharp rise in transonic drag. At Mach number in 

the vicinity of 0.7, the vehicle becomes nwch less stable and remains slightly 
stable at the lowest Mach number. This behavior is a characteristic of 
hypervelocity vehicle and the quantitative results agree with the results 
obtained from a simulation study presented in reference 5. 

The Mach number for neutral point, K = 0, is evaluated from Fig. 9 
at the intersections of the dashed line representing the variation of Cm_- 
and the solid lines for the variations of at different angles-of- 
attack. The evaluation of the Mach numbers which correspond to the stationary 
location of the aerodynamic center can be easily achieved if the aerodynamic 
characteristics are modeled as functions of the Mach number for each angle- 
of-attack of interest. By taking the derivative of Eq.(41) with respect to 
M, we have 


C = C 

ni 


ac, 


Vn 


+ KC, 


M 


^M 


(51) 


When K|^ = 0, substituting back into Eq.(41), we deduce che equation 


Cm - r 

^ "ac . S 
(C„ - C„ )„ ■ CT 


(52) 


which can be solved for two values of M, one near 0.7 and one near 1.2. The 
first value corresponds to a relative minimum level of stability while the 
second value corresponds to the absolute maximum level of stability. 

Hypersonic Trim 

During entry, the trajectory is controlled by the angle-of-attack. 

With only the aerodynamic moment involved, longitudinal trim at = 0 is 
effected through body-flap and elevon deflections, denoted by 6p and 6e 
respectively, with positive value when the trailing edge is down. For the 
exanple CCV, aerodynamic characteristics at hypersonic speed, M = 20.3, 
are presented in Figs. 12-14. Figure 12 plots the normal force coefficient 
Cfi as function of the angle-of-attack for various body-flap and elevon 
deflections. The variations of the moment coefficient as function of the 
angle-of-attack are presented in Fig. 13. To show the trimmable condition 
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Moch number 


Fig. 10. 
function of 


Location of the aerodynamic center as 
the Mach nun^ier (low angle-of-attack) 








as function of the location of the center of gravity, the values of are 
evaluated for a center-of-gravlty at 0.715 of the body length. For any given 
angle-of-attack, with different location of the center-of-gravlty, the new 
Cni curve can be deduced by using the relation (50) written as 


(X2 - Xl) , 

C - c = ( 4 ) C„ (53) 

ffl2 mi ^ ^ c ^ N ' ^ 

where xi/£ and xz/i denote any two locations of the center-of-gravlty In 
percentage of the body length with corresponding moment coefficients and 
Cm . It Is seen In Fig. 13 that at Xcg/£ = 0.715, with ® 

small eleven deflection, negative when a > 42.5® and positive when o < 42.5® 

Is sufficient for hypersonic trim. Based on Eq.(53), If the center-of-gravlty 
Is moved further aft, a positive nose up Increment Is Introduced, propor- 
tional to Cfj and hence large at high angle-of-attack. This will require 
larger positive deflections of both and 6-. On the other hand, a move 
forward of the center-of-gravity can result in a negative pitching moment 
coefficient and a negative eleven deflection Is required for trim. For any 
allowable deflections of the control surfaces. In both positive and negative 
limits, there exists trimmable center-of-gravity range, function of the 
angle-of-attack In the hypersonic range and also function of the Mach number 
in the transonic and supersonic range. 

In general, in the hypersonic range, at any angle-of-attack, for a 
specified location of the center-of-gravlty, the function Cm(a, 5p, 6g) can 
be linearized and we have 




s- + ( 


5 

36 . 


m 


)o 6. 


(54) 


The stability derivatives (3Cm/35p)o and (3Cm/36e)o are evaluated for 
Xp /£ = 0.715 and plotted in Fig. 14 as functions of the angle-of-attack. 
At^any angle-of-attack, for any selected body-flap deflection, 6p, the 
eleven deflection, required for trim is computed from = 0, and 
subsequently the trinmed lift coefficient is deduced from 


\r,. ■ ^ 


3C, 

'<0 ^ jf ) 


(55) 


Similarly, the corresponding normal force coefficient is obtained from 





)o 


(56) 


The stability derivatives in Eqs.(55) and 
positive 5p and 6g for trim increases the 
force coefficient. 


(56) are positive and using 
lift coefficient and the normal 
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For a SSTO vehicle, an aft center-of-gravlty location resulting from 
large engine mass at the rear of the vehicle creates challenging problems 
for the designer. Due to thermal and structural constraints, both the body- 
flap and <?levon deflections have alloMable limits. With these limits, as 
mentl. pf'^. above, there exists trlmmable center-of-gravlty range. To 
asses: '.his range, we follow Freeman and Powell (ref. 5) and consider the 
plot of C«j versus Cn for xrg ■ 0.7151 In Fig. 15. The lowest curve 
corresponds to maximum positive deflections of the control surfaces, with 
■ 10®, 6g - 20®. In general, at any location of the center-of-gravlty, 
for any 6p and 5«, the equation of the Cm versus curve Is obtained by 
eliminating a between the two equations 


'm 


C^(a, 6g) 




(57) 


which results In 



*•111* '•N' 




(58) 


From Fig. 15, In the range of angle-of-attack of Interest, since at maxliman 
control surfaces deflections. Cm < 0, while C,,, > 0 for zero deflection 
angles, the vehicle Is trinmable at this location of the center-of-gravlty. 
For another aft location, from Eq.(53), we have an Increment ACm such that 

AC,„ - CS 


where the constant factor Is given by 





( f ) 


(60) 


with Ax/i expressing the shift In percentage of body length of the center- 
of-gravlty. With Ax > 0, all the ^ curves move upward according to the 
Increment given In Eq.(59). It suffices to analyze the shift of the lowest 
curve for the determination of the most aft of the center-of-gravlty for 
trimmabillty. This Is done In Fig. 16 from which the assessment for the 
example vehicle can be made. 

In general, for any vehicle having similar behaviors for the aerodynamic 
characteristics, consider the Interval of angle-of-attack of interest 
ac[ai, 02 ] which corresponds to the Interval CngCCm. , C^j]. The critical 
center-of-gravlty Is reached when, with maximum deflections of the control 
surfaces, either (^ • 0 for a • ai or when Cm curve has a maximum value 
Cm * 0 at a certain angle-of-attack. 

The first case occurs when 

Cj^(ai) + cC|^(ai) * 0 (61) 

where Cm(ai) denotes the pitching moment coefficient evaluated at aj, with 
maximum control surfaces deflection at reference center-of-gravity location. 
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From Eq.(61) we deduce ; and then the critical location of the center-of*grav1ty. 
Beyond this position, the vehicle Is trlmmable only at high angle>of*attack 
and for each ; beyond this critical value, the lowest value of the trlmmable 
angle-of-attack Is obtained by solving Eq.(61) with ai replaced by a. Of 
course the absolute limit for trlmnablllty Is obtained when this a reaches aa. 


The second case occurs when the following conditions are satisfied 
sitmiltaneously 


Cm(Cj^. «F» *e^ + ■ 0 


K 

c - 0 


(62) 


for a value Cfj within the Interval [Cm , system of equations 

5p and 6 q have their maximum values, and hence the equations can be solved 
for c and Cn* We then deduce the critical location of the center-of-grav1ty, 
through ;, and the critical angle-of*attack, through Cfij. Beyond this position, 
the vehicle becomes trlmmable either at high angle-of-attack or at low angle-of- 
attack . For each ;, beyond the critical value, the two limiting angles-of- 
attack are obtained by solving Eq.(61) with replaced by a. 

Graphically, since for any percentage variation In the center-of-gravlty 
resulting In a value Ci the Increment aC„, from maxlimim deflection Is given 
by Eq.(59), we can plot this equation as a straight line, with a change In 
sign for ; and measure the new value of Cm for maximum deflection from this 
line (Fig. 17). Figure 17-a for a hypothetical vehicle corresponds to the 
first case. At Xcg ■ 0.735t, the vehicle becomes just trlmmable at a * 20*. 

As the center-of-gravlty imives further aft. It Is trininable only at higher 
angle-of-attack. Figure 17-b corresponds to the case of the present example 
CCV. At Xcg ■ 0.739i, the vehicle becomes just trlmmable at Cm ■ 1.350, 
a * 33.1®. ^As the center-of-gravity moves further aft, the vehicle becomes 
trlmmable either at high angle-of-attack orat low angle-of-attack. 
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IV. PHUGOID OSCILLATIONS 


With the static stability assessed, we now Investigate the dynamic 
stability In longitudinal motion. 

In low speed flight, when perturbed from a reference flight condition, 
longitudinal dynamics Is characterized by a phugold oscillation with long 
period and an angle-of- attack oscillation with short period. The phugold 
mode Is a trajectory mode along which the angle-of-atttack remains nearly 
constant while the center-of-gravity of the vehicle oscillates about a 
reference flight path through an exchange between the potential energy and 
the kinetic energy. By extension, we shall use this decoupling to analyze 
the first order effect of density gradient and centrifugal force on phugold 
motion at nearly constant angle-of-attack. 

The phugold, or trajectory oscillations for the three typical cases of 
cruising flight, ballistic entry and glide entry are the subjects of 
Investigation In this chapter. In Chapter V we shall analyze the angle-of^ 
attack oscillations during entry for these cases. 

Cruising Flight 

The simplest motion can be envisaged when the vehicle Is perturbed from 
a steady-state condition. From Eg. (25), at very high altitude, this can 
be achieved by using a small thrusting force to balance the drag and maintain 
a constant orbital altitude. With the thrust aligned with the velocity and 
considering a cruising flight at constant angle-of-attack ao resulting in a 
lift and drag coefficients, Ci and Cn , we have the equations of motion 


m ^ = T - % pSCjj - mg slny (63) 

mV ^ H ~ ^ )cos-y 

At the reference constant-altitude cruise 


R = Ro, V = Vo, Y = 0 
T = Js PoSCp^ V| 

h PoSCj^^Vi = m(go - ^ ) 

It is convenient to use a time transformation 
T '* i V(t)dt 

>'o 0 


(64) 


(65) 
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to rewrite the equations. We notice that this new Independent variable 
denotes the dimensionless arc length. Then, with the prime denoting the 
derivative with respect to t, the system of equations (63), with the constant 
thrust replaced by Its expression In Eqs.(64), becomes 

i 

6C. - ( aSi )sinr (66) 

» Uo ^2 Wo y 2 

y' = iC. - { aBa . ^ )cosY 
Lo y 2 K 

where we have defined the dimensionless atmospheric density 



From the reference steady-state orbital cruise, we have the small perturba- 
tions 

R = Ro + aR 

V = Vo + aV (68) 

Y * 0 + Y 

The gravity acceleration g(R) and the density 6(R) are functions of R. For 
an inverse-square force field we have the linearized acceleration of the 
gravity 


g = g„n - 2 I; ) 


(69) 


The density function 6 Is more complex. In general, by Taylor's series 
expansion 


6 




+ Ro( 


d5 \ ^ 

dR '• R. 


(70) 


where (d6/dR)o Is the first order density gradient. The system (66) can be 
linearized and using the reference condition (64) for simplification, we 
have the equations for the perturbations AR, aV and y 


aR' 


M' -- r R f ^ 1 ^ r M 

U ^Pl •'oV AO fo D “ y 


dR Ro 


Vi 


y' • [c, R,( ^ )„ + - 1 ] |B + 

y2 \i2 


Vi 


4V 

Vo 


( 71 ) 
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Any merit of an analytical study results from the fact that certain parameters 
governing the physical phenomena can be displayed explicitly. From the 
linearized system (71), it is seen that the characteristic roots depend on 
the constant parameters 



p oSRoCis 

Ha 

2m 


a 


2 



(72) 


= la‘( ^ ) * Rfl. ( M ) 
o ^ dR (So ^ dR ^ 


The first term is the drag parameter at the reference altitude. It is function 
of the angle-of-attack selected and the orbital altitude. The parameter a 
is the ratio of the reference cruising speed V® to the orbital speed in the 
vacuum /gofto at the reference altitude. We notice that a -*• 1 as p ^ 0 at 
very high altitude. The parameter a is termed the ^t.ofpheric density 
gradient. From the reference condition (64) we deduce the relation 

«.C. = lijL-iii (73) 

L, ^2 


Hence, ultimately, phugoid oscillation, for any vehicle considered, depends 
on the flight altitude, which provides the values a and 5o, and the trimmed 
angle-of-attack a©» which gives the values and Cl^. The speed parameter 
a is deduced from Eo.(73). The effect of the altitude is entered through 
both the values po (or 6o) and the density gradient a which is net negligible 
at orbital flight. The effect of the orbital speed is due to the fact that 
Eqs.(73) account for the effect of the curvature of the Earth. 

Using the vector perturbation 




we can write the linear system (71) in matrix form 


dr 


AX 


with the constant square matrix A being 


r 


0 


0 


A = 


-(SoCp^a 


-26 oC 


Do 


iLzAll g . J- , 1 i- 
a^ a^ a' 


1 



0 


(74) 


(76) 


(76) 


35 


The resulting stability characteristic equation Is a cubic equation which 
can be obtained by expanding the equation 

det. [A - XI] - 0 (77) 

We have 


1 26 oC- 

X® + 26oC. X^ + -5- C(l-a2)(2-62o)+4*»3x ^ [2-4^-6"a] - 0 

For reason which will be clear later, we define 

^2 , C(l-6")(2-4"g) -K 4*^ 1 
4** 


■2 3 [2-4^-6^a! 


(78) 


(79) 


Numerical values for the density gradient, a, as function of the 
altitude are given In Appendix A. In the altitude range of Interest, a 1s 
oscillating In the vicinity of the value -900. On the other hand, the speed 
parameter 4 varies In the range 0 £4 £ 1, close to zero at low altitude and 
tending to unity In orbital flight. Hence, the positiveness of the parameters 
0 )^ and Is assessed. With the definition (79), we rewrite the 
characteristic equation as 

X^ + 25oC-, X^ + - 26oCn = 0 (80) 

Uo Wo 

This stability cubic equation for the trajectory mode has three roots, one 
real and small root and one pair of complex conjugate roots. The real root 
is obtained by observing that the last term of the cubic equation Is small. 
Hence, we have approximately 


Spiral ° ^*'‘'0, 



(81) 


Since this real root Is positive, hence unstable, the name of spiral mode 
is justified. The exact root, to any order of accuracy. Is obtained through 
series expansion by putting 


^spiral “ ^ I 

where a is a constant near unity to be determined. 
Eq.(80) we obtain 


(82) 

By substituting into 


25 C 

a = 1 - ( a“ [1 M ^ )^a] (83) 

U) 0) (i) 

The equation is In the form 

a ’ p + ef(a) (84) 
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where, with a sinaU c, Lagrange expansion 


* ■ p ' i >"■’ '®=> 

is dictated. Here, the small term (26oCQ^/tu)^ plays the role of e. Hence, 
we have explicitly 

25 C 

+ ( — ^ )' ( ^ )‘ [1 + ( ^ ) n [2 + 3 ( ^ )^] ( 86 ) 

0) 0) 0) (1) 

- { ( M‘ [1 + ( § )^][l + 2 ( ^ )n[5+ 6 ( + ... 

U) U) Q) (U ui 

Since a = 1 already gives a good approximation, this explicit series solution 
Is very accurate. Next, by factorizing the cubic equation, Eq.(80), using 
the real root (82), we have the equation for the remaining pair of complex 
conjugate roots 

+ 26oCn [1 + ( ^ )^a]x + j = 0 (87) 

Uq (i) a 

where, of course, the value of a is given by Eq.(86). From this equation we 
have for the phugoid mode which Is a damped oscillation 


and 


•"^^phugold' " “ 




)“ [1 + ( 5 

li) 


Using Eq.(83), we have another form for the phugoid frequency 


(89) 


>”(Shugoid> * “ V’ - < Cl t ( I )=.][! - 3( I )^a] (90) 

Returning to the dimensionless time, t, as defined In Eq.(65), If the small 
speed perturbation Is neglected. It 1s simply 


T 



t 


Then, the time to double for the spiral mode Is 


(91) 
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(92) 


^double poSCplpa ^ 5 ^ 

The spiral mode is a new mode Inherent In near orbital flight. Since we 
assume a small constant thrust which exactly balances the drag In cruising 
flight at constant altitude, any perturbation causing an Increase In the 
altitude with smaller drag will create an excess of thrust which has the 
effect of increasing the radial distance, hence sending the vehicle Into an 
outward spiral. The reverse effect Is true for any decrease in the altitude. 

Now, since the phugold mode Is a damped oscillation, we consider the 
spiral effect on the change in the flight path angle and the speed. We 
start with the solution for the radial distance 


Ro 


ee 


Xt 


(93) 


where x is the spiral root and e Is the initial perturbation with e > 0 for 
an outward spiral. From the first equation In the linear system (71), we 
have 


y = eXe^^ (94) 

which shows that y varies as 4R. On the other hand, from the last equation 
of the system, we have 

^ = f - 0-4^)c + 4" - (95) 

Since x^ is small and o -900 (see Appendix A) while a 1 when p -► 0, the 
bracketed coefficient is negative at very high altitude and is positive at 
low altitude. Hence, at very high altitude, the speed varies In the opposite 
direction as the radial distance, that is the speed decreases along an out- 
ward spiral and Increases along an inward spiral. At lower altitude, the 
speed varies In the same direction as the radial distance. Neglecting the 
small term x^4^, the condition for the speed to vary In the same direction 
as the radial distance Is 

-d-4^)0 + . 2 > 0 (96) 

In terms of as defined In Eq.(79), we rewrite this condition as 

(4)24** > 2 (97) 

In a more explicit form, using Eq.(73) to calculate the speed parameter 
we rewrite the condition (96) 

-{2+o)6oC, > 1 

Lo 

since -o Is large, we simply have 
Lo 


38 



If the equality sign is used in this equation, upon solving we have the 
altitude where the speed Inversion occurs. Graphically, since the right- 
hand side of the condition (98) Is purely a function of the altitude, for 
any planetary atmosphere, the function can be plotted versus the altitude. 
Then, for any wing loading condition (W/S), and lift coefficient, C|_^, the 
altitude for speed Inversion for the corresponding vehicle can be assessed. 

Concerning the phugold oscillation, from Eq.(88), we have for the 
damping. In real time 


_ 2m log2 1 


( 99 ) 


The period of oscillation Is 


6oC, 


P . ^ (1 . ( - 1 ^ )“ [1 + ( f )"a] [1 - 3( ^ 


b> 


Expl lei ting for u, we obtain 


( 100 ) 




( 101 ) 


where Lo Is the lift force in reference cruise. Then, the period of 
oscillation In phugold motion has the final form 


P = 






6oC 


(1 - (-^)ni + (^)'a][l - 


(102) 


If we neglect the small term (SoCp )^i and take the second bracketed term as 
unity we obtain an expression that*1s Identical to Laltone and Chou's 
(Ref. 12). Eq.(102) which takes Into account the effect of the density 

gradient and the centrifugal force correctly predicts the phugold period in 
the entire range of cruise altitude from low altitude for airplane dynamics 
to high altitude orbital flight. Asymptotically, in orbital flight with 
So -► 0, Lo 0, Vo goRo the period tends to 


P (orbit) = = 27 t /RiTgl (103) 

9o 

which is precisely the orbital period along a circular orbit at distance Ro. 

The relative balancing effect between the speed and the density 
gradient on phugold period can be displayed by considering the simplified 
formula 
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p 





(104) 


For flight at very high altitude In the near vacuum, vS/goR© Is near unity 
and, although -a Is large, because Lo/mgo Is near zero, the ef'ect of 
orbital speed Is predominant. At medium altitude, as V© decreases and L© 
becomes larger, the effect of the density gradient becomes more pronounced 
and should be Included. At lower altitude, for airplane-type dynamics, we 
write the bracketed term 


^ ^ ^ m^ ^g©R©^^^g©R©^ “ ^mg©'^'^^ 


(105) 


At sonic speed, taking V© « 340 m/sec, /g©Ro * 7919 m/sec, we have V|/g©R« * 
1.843 X 10”^ while Lo/mg© Is near unity. Hence, for all practical purpose, 
we can take 


{ 




(106) 


Then, by setting L© = mg© as condition for equilibrium steady-state flight, 
a condition which can be deduced from the last equation In system (64) by 
using the same approximation, that Is with V|/g©R© << 1, Eq.(104) becomes 


P (airplane) 



or 


P (airplane) 



2ir 


2qo 1_ ^ X 

y2 ■ P dH ^ 


(107) 


It Is shown in Appendix A that the density gradient term Is equal to 
-1.577 X lO'^^/m for an Isothermal atmosphere while for sonic speed, M = 1, 

V© * 295 m/sec In the stratosphere we have 2g©/vS » 2.255 x 10"‘^/m. The 
speed and density gradient effects are now of the same order of magnitude. 
For very low speed, the density gradient effect can be neglected and we have 
the classical formula 

P (airplane) = (108) 

9o 

which, as well known, gives the phugoid period within 10% of the value 
obtained by considering the full set of linearized equations taking Into 
account the coupling between the phugoid and the angle-of-attack modes. 

From the physical discussion above, we have seen that the density 
gradient effect Is most pronounced at the Intermediate range of speed. 
Mathematically, since 
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mo" " ^ • 1 - 

IHQo 9o"o 

we express the bracketed term 1n Eq.(104) as 
{ /a** + (l-a*)(2-4^a) 

Vie have seen that near orbital flight, we have asymptotically 

{ « wi » a* (m ) 

while at very low speed we take 

{ “ /2 (112) 

The thr ee fun ctions u above are plotted versus the speed parameter 
& • Vo/^goRo In Fig. 18 using an average value -o • 900. It Is seen that 
the Influence of the density gradient Is Important and has the effect of 
shortening the phugold period. The approximation wi and B 2 can only be used 
near the limiting cases either near orbital flight or at very low speed. 

Ballistic Entry 

We now consider the Important case of ballistic entry. With T ■ 0 and 
Cl “0, the system of equations (66) Is reduced to 

f . -«C. - (113) 

V Up y2 

y' a - -§^)C0SY 

y2 R 

In the second equation of this system, the drag force during entry Is large 
as compared to the small component of the gravity along the tangent to the 
flight path. Furthermore, we shall use the approximation Ro/R = 1 In the 
last equation. Then, we have the simplified system 

R* 

r7 “ 


(109) 

( 110 ) 




(114) 


It has been verified in reference 13 that the numerical solution of this 
system Is accurate as compared to the solution obtained from the Integration 
of the exact equations. Following the analysis In that reference, we use 
the Chapnan-Yaroshevskll variables 
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(115) 


pSCn 

Y - 

^ - /BRo slny 

X « log (gRo/V^) 


where 6 Is the Inverse of the scale height In a locally exponential atmosphere, 
that Is 

dp ■ -"^pdR (116) 

Then, we have the equation for the altitude which 1s here represented by the 
new dimensionless density Y 

Y' - -iRoYsIny (117) 

For the flight path angle which Is represented by the variable "I*, we have 

--/fRi (1 - -3^)cos^ (118) 


Changing the Independent variable to V and then to x, we have successively 


dV 



dV 


W <’ - 


and 


dx 




. (e^ - l)cos^Y 
dx Y 


(119) 


( 120 ) 


Because of the physical constraint on the deceleration during entry, ballistic 
entry Is generally effected at small flight path angles. Hence, we have the 
dimensionless equations for ballistic entry at small and medium angles by 
, that Is 


( 121 ) 

S “ Y ' 


taking cos'y - 1 


dx ^ 




t. 


For ballistic entry from circular speed, the Independent variable x Increases 
from the Initial value zero to about 6.3 when the speed has decreased to the 
order of the speed of sound at low altitude. On ti^ other hand, the altitude 
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variable Y, which Is proportional to the density, starts with near zero value 
at the top of the sensible atmosphere. For iwst vehicle physical parameters, 
the value of Y at sea level is near 60. Expressed In the change In the 
altitude, we can consider the change In an exponential atmosphere 

IaH-1o9 7 (122) 


where Yi Is the Initial value of Y and AH * - H Is the variation in the 

altitude from the initial height. This formulation is convenient In that It 
Is Independent of the physical characteristics of the vehicle. 


Accurate series solutions for the system (121) have been given In 
reference 13. For the case of ballistic entry from a decaying orbit such as 
the case of the reentry of the skylab, we have the initial condition 


v( 0 ) • 0 . 

0(0) » 0 


(123) 

with the solution 




V - -2- 
/3 

[1 1 ^(|) . . ^(|)> 

. 20021 

(124) 

and 




0 ■ /3 [1 

* |(t) * ^ 

20021 ,x,4t 
165240"^?^ J 

(125) 


For entry with non zero initial flight path angle, we have the initial 
condition 

Y(0) ■ 0, 0(0) * - /Wo slny^ » c (126) 


where Is the initial flight path angle. From the definition of a in 

Eq.{72), with a locally exponential atmosphere as defined in Eq.(116), It is 
seen tjiat "^Ro * -a. Hence, as shown in Appendix A, we take the average 
value BRo • 900 for evaluating the Initial parameter c. 

For the solution In this case, we write the system (121) in the form of 
Yaroshevskii non-linear equation 


Y 



X + 


X' 



♦ 


* • • 


027) 


and evaluate the formal solution 


Y 



satisfying the Initial condition (126). It Is found that 


(128) 
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*. • c. ^ . *. ■ (1 - • 7K - -7 * "T > 

' ' ' (12») 


To the same order, we have 

09 

^ • I b.x*^ 
k-0 


(130) 


where 


h« • c, bi ■ ^ . b 2 


Jc 




c‘ c' 


b% 


10 ^ 20 17 


96c 


(1 . -L«. + ^ ) 

3c^ 3c'* 3c‘ 


(131) 


The variations of the flight path angle for entry at circular speed 
with different Initial flight path angles are presented In Fig. 19. The solid 
lines Indicate the exact mxnerical solutions of Eqs.(121), while the dashed 
lines represent the present analytical solutions. Solution for entry from a 
decaying orbit Is very accurate. On the other hand, for non-zero Initial 
angle, the solution Is only accurate when > 1.5*. For very small Initial 
angles providing small values of c. because of small radius of convergence 
of the series solution, the solution Is restricted to the high speed range. 

We now Investigate the small perturbations from a reference solution 
Yo. ♦o. Let 


Y - Yo + AY 


By substituting Into system (121) and sln^Ufylng, we have 





dM . («* - Day 

dx ■ V,(Y. ♦ 4Y) 


032) 


033) 


Hence, the perturbation in altitude Is obtained from the non-linear 
differential equation of the second order 


dfw . te* - Day 
Y.IY, .AY) 


• 0 


034 ) 


where the function Yo(x) used Is the reference solution, t*'st Is the solution 
(124) for entry from decaying orbit and solution (128) for entry with non- 
zero Initial angles. By the change of variable 


(135) 


y 


/x 
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Fig. 19. Variation of the flight path angle 
for ballistic entry at circular speed. 
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0 


(136) 


we obtain the new non-linear equation 

dx^ ^ 4x^ Y|(l +y/x/Yo) 

Going back to the initial system (121), it Is seen that for large x, 
eX _ 1 = and the asymptotic solution for Y Is 

Y = 2e^''2 (137) 

If this solution Is used for Yo In Eq.(136), we have, with the same 
asymptotic approximation 

+ + 1 ^=0 (138) 

dx^ * 4x^ 4(1 +y/x/2a*/^) 

Hence, It appears that y Is a damped cscillatlon, and In the linearized form 
we have 

dx“ X >)x 4 ^2 

This is a special case of Bessel equation and the general solution Is 
obviously 

y = -L [Cicos ^ + Czsin 4] (1^0) 

/x ^ ^ 

where Ci and C 2 are two constants of Integration. In general, since y Is a 
damped oscillation, Eq.(136) can be linearized and we have 


. J-]y.O (141) 

dx" * v; 4xX 

For phugold motion, perturbed from entry from a decaying orbit, we use 
Yo as given by Eq.(124) with three terms 


V 

• O 




(142) 


Substituting Into Eq.(141) and expanding the coefficient of y for small x, 
we have 


£/ + 1 ^ . J_ (2 . X 

dx‘ * 4x“ 




(143) 


It can be shown that this equation can be transformed Into a confluent hyper 
geometric equation (ref. 14). If the term x^/4 Is neglected, we have the 
equation 
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( 144 ) 


+ 1 +l2ixli . Q 

dx^ * <'* 4x» 


which can be transformed into a Bessel equation of imaginary order (ref. 15) 
For very small x, if we consider the coefficient of y as simply l/2x^, we 
have obviously an Euler equation with general solution 


/ /2 X / /2 X 

y = CiCos( log X ) + Czsin ( *2- log x ) 


(145) 


where Ci and Ca are constants of integration. It is enlightening to derive 
this solution by considering the general equation (144). By the transforma- 
tion of variable 


X = 


the eouation becomes 




.2 C dc 


(146) 


(147) 


d5^ ’ ’’ 

with general solution 

y = Cl J^/2 ^ 

where the functions Op(x) are Bessel functions of order n and argument x. 

2 ^\\ ' 2(2n+2) 2-4(2n+2)(2n+4) " 

To obtain the solution in the real form, we consider the two linearly 
independent solutions, putting n = ±1/2 , and changing x into/x in the 
expression (149) 




4(l+1/2) 32(1+1/2 )(2+i/2 ) 

v3 

X ^ + . ..] 

384(1+1/2 )(2+i/2 )(3+1/2 ) 


4(l-1/2 ) 32 0-1/2 )(2-1/2 ) 

y3 

2 + . ..] 

384(1-1/2 )(2-i/2 )(3-i/2 ) 

After some manipulations 
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(150) 


y. o + 1g.(x)] 

yx ^ x-'’^2/2[f.(x) - 1g,(x)] 

where 


fi(x) * 1 + 12672 ‘ 

„ f^\ /2 X X . x^ \ 

9 i(x) = (1 - 16 ^ m ‘ 


(151) 


Since a linear combinction of solutions is a solution, we combine linearly 
the solutions (150) to obtain new solutions in the real form 

^~2 /Z 

yi = fi(x)cos( -j- log x) - gi(x)s1n( -j- log x) 

(152) 

/2 /2 

yz ’ fi(x)sin( Y ^09 x) + gi(x)cos( log x) 

The solution (145) is obtained using the approximation fi(x) = 1, gi(x) = 0 
for small x. Using the two linearly independent solutions (152), with the 
definition (135), we obtain the solution for the perturbation AY in the form 

AY = A /x [fi(x)cos( Y ^09 - 9 i(x)sin( y ^°9 jjl (153) 

where A and Xo are now arbitrary constants of integration. By the first of 
the equation (133), we obtain the perturbation A$ in the form 

A$ = -At [f 2 (x)cos(-^ log A " 92 (x)sin(-~ log (^54) 

t Ao t Ao 

where 

92(x) = /I (1 + I - 2^5344^ 

The solutions (153) and (154) with the functions fj(x) and g^(x), i *1,2, 
not only provide the solution for the phugoid oscillations for entry from a 
decaying orbit but at the same time yield accurate analytical solutions for 
entry at very small angles. We recall that, as stated above, the existing 
solutions (128) and (130) cannot be used when the value of c is small. 

It is now seen that, with a misalignment in the entry angle, yi t 0i 
the difference in the altitude of the perturbed motion, as compared to 
entry at yi = 0» increases. This is to be expected since entry at non -zero 
initial angle has the effect of drastically decreasing the range, providing 
a large difference in the altitude. On the other hand, the variation in the 
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flight path angle, as shown explicitly in Eq.(154), is a danced oscillation. 
Hence, in this sense, phugold oscillation for near grazing ballistic entry 
is stable. To show this quantitatively, we use Eq.(154) to compute the 
perturbation from the reference solution as given by Eq.(125). The constants 
of integration A and Xo in Eq.(154) are evaluated using the initial condition 
(126). Because of the singularity at x * 0, the constants are evaluated at 
a certain point beyond the entry point. By convention, we take this point at 
X * Xo. For the computation at this arbitrary initial point, we refer to 
the definition (132) and use the following approximations; 

Since in this case c 0, the functions Y and $ are initially given by 
the solutions (128) and (130) which are valid for small x, even in the case 
of small c. The functions Yo and $o are qiven by Eos. (124) and (125) and 
the functions AY and A4> are given by Eqs.(153) and (154). Then at the point 
X = Xo, and to the order of x*, we have 


^ . |i) + A/x fi(x) = cx + 0 - ^)x’ 056) 


/3 


and 


5x . 7x 


1 . 1 


1 


36 288^ 2/- ^ c^ ‘‘‘ 4c^^ ’ ^ 2 ^^ 

where, from the expressions for fi(x) and f 2 (x), it suffices to take 
fl(x) = 1 - 
f 2 (x) = 1 - -fl 


(157) 


(158) 


In these equations, subscript zero in Xo has been omitted for convenience. 
Then, by solving the system of equations (156) and (157), the initial value 
Xo is obtained from 


(4 

/3 


+ X + 


T 


288 


) = c(l + 


4^ 


^(3 

2c^'^ 


12 ' 





and subsequently the constant A is given by 


(159) 


/3 ® 32 ^ 


_1 

/X 


[c + 


3x . 5x' 


2c 12c 


(1 - — )] 


(160) 


where we recall that, in the two equations above, c = - /ERo siny^ * 

-30 sinvi*, and x = Xo- The solution for A$, computed for -yj in the 
interval from 0® to 2®, displays, as expected, a damped oscillation with 
long period. Furthermore, when added to the nominal value <I>o{x), it gives 
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accurate solution for shallow ballistic entry. The result is plotted in 
Fig. 19 and is indistinguishable from the exact numerical solution. 

Half of the period of oscillation can be defined as the Interval 


between two consecutive zeros, xi and X 2 of A$, as given explicitly in 
Eq.(154). Let 

tan[i^(x)] » g 2 (x)/f 2 (x) (161) 

Then, the zeros of are also the zeros of the function 

A$ = cos[^ log ~ + \i>(x)] (162) 

Hence, if Xi is the first zero and xz is the next one, we have 

^ log ^ + i|^(x 2 ) - i^(xi) * TT (163) 

Similarly, we define 

tan[iF(x)] = (164) 

For the oscillation in the altitude, if xi and xz are the first and the 
second zero of the function AY, respectively, we have the relation 

log + ^(xz) - = TT (165) 


Figure 20 plots xz versus xi for the Eqs.(163) and (165). The value xi 
increases from zero as -Yi Increases from zero. Since the practical range 
of X is between 0 and 6, except for near zero values of xi, which correspond 
to very small values of Yj , the value o/ Xz is large. Hence practically the 
plots of the variations of the flight path angle for Yi ^ 0 cross the plot 
for the case of orbit decay, y^- = 0* only once (Fig. 19). 

Next, we consider the phugoid motion as perturbed from entry at medium 
and large angles. As seen in Fig. 19, with an initial perturbation 


A(}>(0) = - /BRo (sinY^ - sinY^Q) = Ac (166) 

where y-|q is the initial reference angle and Y-f = Yio is the initial 

perturbed angle, the function - Ay, and hence A$, is a slowly decreasing 
function. Again, in the (-y. V) space, the two neighboring trajectories 
cross each other once at low speed. Hence, the phugoid motion is a slowly 
damped oscillation with long period. Since accurate analytical solution has 
been given in Eqs.(130) and (131), it suffices to take the difference 
between two neighboring solutions to obtain the solution for phugoid motion. 

To display the damping characteristic of the phugoid motion, we consider 
the linearized perturbed equation (141), with the reference solution Yo(x) 
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Fig. 20. Relation between the first two zeros 
of the functions AY(x) and A<I>(x). 




given by two terms of the series (128). Then we have the equation, to the 
same order as Eq.(143) 


dx" 


+ 1^ 
X dx 


2c- 



-L]y = 0 
4x^ 


(167) 


Again, this equation can be transformed Into a confluent hypergeometric 
equation. For medium and large Initial entry angle, c^ Is relatively large 
and we can neglect the constant term, which Is small. In the coefficient of 
y. Then, by the transformation of the independent variable 


C = 


2 /x 
c 


(168) 


we obtain a Bessel equation of order 1 


. 1 ^ + (1 


— )y - 0 


(169) 


From the definition (135), we have the solution 
4Y [CiO.(^) + C2V,(^)] 


(170) 


where the general Bessel function J^(x) Is given by the series (149) and 
Yi(x) is Neumann's Bessel function of the second kind of order 1. By 
considering the relation 

Yi = Jilogx--^“-jJi+-^J3-... (171) 


and from the condition AY(0) =0, It is clear that C 2 * 0. From the first 
of the equations (133) we have the solution for A$. Using the relation 

CJI(C) * - Ji(c) + ?Jo(;) (172) 

for the derivative of the Bessel function with respect to Its argument, we 
have 


A* - J. (173) 

Using the Initial condition (166), we have the final solution for the 
perturbation in the flight path angle 

A4 > AC 0. (174) 

V 

since Bessel function of order zero is a damped oscillation, the charac- 
teristic behavior of phugoid motion perturbed from ballistic entry Is 
established. Furthermore, Eq.(174) gives accurate solution for 4>(x) when 
added to the reference solution $o(x). Since the first zero of Jo(^) Is 
i;i » 2.40482... two neighboring trajectories. In the (-y, V) space cross 
each other for the first time at the value 



X| - 1.4458 c2 


(175) 


For medium and large entry angles, this corresponds to a very low speed. 
Glide Entry 

For glide entry, we consider the system (66) with T « 0. By neglecting 
the small component of the gravity alon^^ the tangent to the flight path and 
using small angle approximation we have 



f- • - (176) 

y' • 6C, - - 1) 

y2 


In the last equation, the approximation Ro/R = 1 has been made. Then, using 
the density function Y as defined In Eq.(115) and the new speed variable 


V * 


Jd 

gRo 


(177) 


as the Independent variable, we obtain the equations for glide entry at 
small flight path angle 


<iy - Y 

dv V 



--(-) + 
2 ^ 0 ^ 


1 


/^Ro Y 



(178) 


For small flight path angle, the dimensionless arc length t, as defined In 
Eq.(65), Is the same as the longitudinal range. From the second equation In 
system (176), we have 


dx _ 1 

/IrI Y V 


(179) 


For the value of the llft-to-drag ratio, L/D » 1. 1 
Integrated for glide entry from circular speed, V > 


5, the equations are 
1.0, and an entry altitude 


which corresponds to Y^ ■ 0.001 with an Initial angle y* * 
tion In the altitude during the glide Is shown In Fig. 21. 
displays the oscillatory character of phugold motion. For 
of the altitude, we can use the equilibrium solution which 
fact that the variation In the flight path angle Is small. 
Hence, from the second equation In system (178), we have 


-4®. The varla- 
The figure 
an average value 
expresses the 
that Is dy/dv ^ 0. 


54 



This average solution Is plotted In dashed line In Fig. 21. To have an 
explicit solution of the deviation from this line, we use the new dependent 
variable 


2 . ( !»/.Piyi - 1 (181) 

2(l-v) 

As a result of this definition, except at the beginning of the entry 
trajectory where 2 can be large, the perturbation In z Is a small quantity. 

By substituting Into the system (178), we have the equations 

dz _ -gR» (L/D)Y , (z-H) 

dV 2(1-7) v(l-7) 

(182) 

, -Ik/Jll 

dv 2v(z+l) 

By eliminating y between the two equations, we have a second order non-linear 
differential equation for z 

7(1.7)^ . (1+7)^ + « 0 (183) 

dv^ dv 4(z+l) V 

The equilibrium solution (180) does not provide the oscillation in the 
altitude but gives an average value with good accuracy. Hence, the function 
z gives this correct oscillation and tends to zero near the end of the 
trajectory. By linearizing Eq.(183) we obtain 

7(1-7)^ - (1+7)^ + [BsJL/Oli + ±32 - - ^ (184) 

dv^ dv V V 

Since "HRo - 900 Is large, In the homogeneous equation, when v Is not too 
small, we can neglect the term 1/7 and obtain a hypergeometric equation. 

In general, with the change of variables 

2 y , V » cos^U , 0 < y < 5 (185) 

(tany)3/2 

the linear equation (184) Is transformed Into 

. CER. iky * ^ ^ ]U ■ - * {,8j, 

dy^ ° 4 sln^y cos^y cos*y 

In the homogeneous equation, the non-constant term In the coefficient of U 
Is 


I 


J - < . -Ms— (187) 

4 sln^u cos*vi 4v(1-v) 

In the range of speed of Interest, when v decreases from 0.95 to 0.01, this 
coefficient Increases from a negative value -14.737 to a positive value 
24.242. On the other hand, the term FRo(L/D)^ Is very large. Hence, the 
solution of the homogeneous equation In z Is practically 

z • — 3 {C, cos|/Br7 (4)>i] + C2 !ln&^ (t)u]> (188) 

(tanu)S/2 '' 

To this equation, we add a particular solution of Eq.(186). Considering the 
solution (188), when v -*• 0, tanu -► » and z -► 0 as expected. Furthermore, z, 
and hence the altitude Y by Eq.(181), has an oscillatory motion In y with 
frequency /"PRo (L/D). Actually, this should be viewed as an oscillation 
about the average altitude as given by the equilibrium relation (180) which 
has been used as the reference solution In the linearization. When y varies 
from 0 to tt/ 2, th e ar gument of the trigonometric functions In Eq.(188) 
varies from 0 to /IPRo 7 t(L/D)/ 2. Hence, the number of oscillations Is 
approximately 
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L/D « I 



Oscillations In altitude during glide entry. 






V. ANGLE OF ATTACK OSCILLATIONS 


For the angle-of-attack oscillation during entry, we consider the full 
set of equations, rewritten here for convenience 

II ■ V slny 

dv T 

■ S ■ -TUT ■ « 

dv “5C.V' 

* It • -TT" - <9 - X )cosT 

do . 3g , ’x ■ ’t , , 

3t ■ ~1T^ ' ?R < ~T~ > 

de . . . V 
3t ■ 9 * l( 


(190) 


0 • Y + a 

The elimination of 0 and q results In the following exact equation for 
the angle-of-attack a 


dt^ 


SC. V 


dC, pScC„V* 


I -I. 


* -?r H * ^ - -TT- ^ ‘ ^ 


PSC, 


* - !f>9 slny cosy - g cosy - (f|)" W 


(191) 


pSC 


* i (-5=^ + cosy) - 0 


m ‘ zm 


In deriving this equation we have used an Inverse-square gravitational 
attraction to preserve the correct behavior of the motion In orbital flight. 
If the time transformation (65) Is applied, we have 


T 


1^ i 




( ) 





( )•• * (^)( ) 
Ri 


(192) 


where the prime denotes derivative with respect to t. Then, upon using In 
Eq.(191) we have the non-linear equation for a 
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a 


N 


+ 


11 

V 


a' 


+ 




)sin2(Y+o) + - ^)(-SS*)s1nY 


cosy 


09i) 



RopS C| 


+ -SfiA cosy) 
V* 


0 


By the display of this equation, It Is natural to use the following 

I 

k (194) 

Notice that the density 6 and the kinetic energy v are the same as have 
been defined In the equations (67) and (177), respectively. Then, we have 
the equation In the final form 

_ C 

^ a’ ♦ 5'C, + 5C; - 6(^) ^ (|a)sin2(y+ot) 

2v *■ '• c y 2v ^ 

6C 

+ - 4-)(:lr)s^ny cosy - ^ cosy - fiH, C« (195) 

. U/jal (JC, . . 0 

V ^ V 


dimensionless variables and parameters 


It Is more enlightening to linearize this equation separately for the 
different cases of reentry considered. 


Cruising Flight 

To show the correct behavior of the angie-of-attack oscillation In the 
asymptotic case of orbital flight, we consider the case of circular orbit 
In the vacuum with T ■ 0, 5 ■ 0, y ■ 0, V ■ 1 . Then the non-linear equation 
(195) is reduced to 


a" + s1n2a » 0 


(196) 


This Is the equation of a compound pendulum, and while In orbit. In Its 
stable equilibrium position, the vehicle points Its axis of minimum moment 
of Inertia toward the center of the Earth. The equation Is Integrable by 
the use of elliptic Integral. Since In this case, x ■ V.t/R*, the period of 
oscillation. In real time. Is 


where ao is the Initial, perturbed angle-of-attack. Neglecting order of al 
and higher, and noticing that Vo “ /goR© » upon comparing with the expression 
for the phugoid orbital period as given In Eq.{103), we see that the ratio 
of the two periods Is 

098 ) 

Depending on the Inertia characteristic jBj, the period for the angle-of- 
attack mode can be greater than the period for the phugoid mode. Hence, 
the names of long and short period modes, traditionally employed In air- 
plane longitudinal dynamic stability analysis, are no longer appropriate in 
reentry dynamics. For the example SSTO vehicle with physical characteris- 
tics given In Appendix B, we have |B| = 2.3596. Hence, the period of 
oscillation for the angle-of-attack is equal to 65% of the orbital period. 


For cruising flight at lower altitude, aerodynamic moment is restored. 
Using a thrust to balance the drag and the lift to maintain constant 
altitude flight, we recall the equilibrium condition (64), now written as 





(199) 


On the other hand, with a perturbed angle-of-attack a, we have the 
linearized expressions 



+ Cd 0. 

a 


C, = C, + C, a 

L Lo 


(200) 


while for the pitching momemt coefficient Cm, experience has shown that we 
must express it as 

® (Y+u) (201 ) 

mm m* m 
a a q 

where q is the vehicle pitch rate referred to the local horizontal. In all 
the cases considered, ^either cruising flight, ballistic entry or equilibrium 
glide, the condition y = 0 is enforced. On the other hand, the stability 
derivatives and Cm_ written in dimensionless form are based on the body 
length, or wing chord^ c as reference length. Hence, from the term in Cjj, 
in Eq.(l95) and using the time transformation (192) for the derivative 
we consider 




3C. 


(ila)c (i^)C a + (^) (^)a' + (^) ^ 


c ^ 


c a 


8 a 




Changing notation, we write 


(^)C = (^)C„ a + C„ a' + C„ a' 

' - ' m T m m* m^ 

c c a a q 


( 202 ) 
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with the definition 




c 9 a 



(203) 


Using the equations above In the basic non-linear equation (195), and 
noticing that 7 = 4^, R = R©, y * 0 and linearizing a, we have the equation 
for the angle-of-attack oscillation In circular orbit 

a" + 2ria’ + n^a = 0 (204) 


where 


n = 


n' = 


- rc 

2 


(C + C ) 

' nr m ' 
a X.] 


a 


T 


_B_ 6(2-^") p f Rfl-N 5 r 

4^ 6^ a c y a 


(205) 


The viscous damping term is always positive. Hence the condition for 
stability of the angle-of-attack mode Is that n^ is strictly positive. By 
Inspection of this term, for slender type vehicle B < 0 such as In the case 
of the present SSTO vehicle. On the other hand, as has been analyzed In 
Chapter IV, for cruising flight the speed parameters Is In the range [0, 1], 
with near zero value for airplane flight mode at low altitude and close to 
unity near orbital flight. Hence, the negative effect of the second term 
Is diminished with smaller drag curve slope Finally, the stabilizing 

effect is due to the static stability derivative Cma which must be negative 
for stability in cruising flight. Otherwise, thrusters control must be 
Introduced. With a negative value for Cma, ® smaller pitch moment of 
Inertia, ky, contributes more to stability. It Is obvious that a large 
value of 6 contributes to the positiveness of the term n^. Physically, this 
means that, at low altitude a restoring aerodynamic moment overcomes the 
pitch down effect of the gravity gradient. The altitude, below which the 
vehicle is stable. Is obtained by solving the equation n^ = 0, that Is 


_B — 5(2-a^) C = 0 

4 4 a c y a 


(206) 


This altitude s very high so that we can use the approximation 4=1, and 
also neglect the small term In Cg^ to rewrite the equation as 


Rip 


6(1 -I ) 
X z' 


ScC. 


a 


(207) 


which separates vehicle and atmospheric characteristics. For the vehicle to 
be stable at nearly all altitudes, that Is for very small value of p In 
this equation, this requires large negative value for Cma small difference 
In moments of inertia about the axes In the plane of syimetry. The left- 
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hand- side of this equation Is plotted versus the altitude In Fig. 22, while 
the right -hand-side Is used to locate the limiting altitude for stability 
for any vehicle. The small circle In the figure represents the present 
example vehicle. 

In the altitude range where the vehicle Is stable, from the damping 
term n and with the definition of t, we have the time to damp to half 

* 4m log 2 

"^half (C^. t ) (?08) 

S.SV, [c^ — 

a y 

It Is seen that small wing loading, low altitude, high lift curve slope, 
strong damping In pitch and small pitch moment of Inertia, all contribute 
to shorten this time constant. 


At high altitude where n is small, the frequency for the angle-of-attack 
oscillation Is practically n. Hence, we have for the period in real time 



(209) 


With the approximation 4^ = V|/goRo = 1 , we have the expression 


P = 


2tt Vq 


9o[B - 6Cp 

a 



(210) 


As a consequence, near the altitude where the vehicle becomes unstable, 
n - 0 and the period for the pitching mode 1s Infinitely large. Since the 
period of oscillation for the phugold mode is finite, with the limit being 
the orbital period, and the period oscillation for the angle-of-attack mode 
is short at low altitude, there exists an altitude where the two periods 
are equal. Cruising at this altitude, the vehicle Is subject to unstable 
resonance effect. The resonance altitude was detected by Etkin (ref. 16) 
and the resonance effect was discussed by Vinh and Dobrzelecki (ref. 17). 

It has been shown In Chapter IV that the phugold frequency in cruising 
flight Is 00 as given by Eq.(79). Hence, the resonance altitude is obtained 
by setting n = oo. This altitude Is high so that we can use the approximation 
4=1. This amounts to simply using oo = 1. Neglecting the term CQci 1n the 
expression for n, and separating vehicle and atmospheric characteristics as 
has been done In the case of limiting stable altitude, we obtain the equation 


2[3(I -I )-I ] 
R^oP = 


Sc C 


m 


( 211 ) 


Plotting the right-hand-side of this equation In Fig. 22, we have the 
resonance altitude. It Is obvious that this altitude Is below the limiting 
altitude for stability. If a strictly exponential atmosphere with height 
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scale 1/1 Is considered, by taking the ratio of the two equations (207) and 
(211), we have the difference between the resonance altitude and the 
limiting altitude for stability 

log (5^) (212) 

This Is quite an Interesting formula since the altitude difference Is purely 
a function of the Inertia property B. For the example SSTO vehicle, we have 
BAH* 0.3533. Using Chapman's approximation for the atmosphere, BRo = 900, 
we can write the equation as 

With an average value Ro = 6467 km at the flight altitude, for the vehicle 
considered, we obtain AH* 2.539 km. 

A final note to close this section Is In order. For the phugold mode, 
as has been discussed In detail in Chapter IV, the period, customarily 
called the long period. Increases from small values, of the order of a few 
minutes, for airplane- type cruise at low altitude and tends asymptotically 
to the orbital period when the atmospheric density tends to zero at very 
high altitude. On the other hand, for the an gle-of- attack mode, the period, 
mistakenly called the short period, also Increases from very small values, 
of the order of a second, to infinity at the unstable altitude where the 
gravity torque and the aerodynamic moment are equal. This altitude is 
above the resonance altitude, where the two periods for phugoid and angle- 
of-attack modes are equal, by an altitude difference as given by Eq.(213). 
The motion is unstable above the limiting altitude but after a pitch down 
to point the x-axis along the vertical, the gravity gradient torque becomes 
predominant and the pitch motion is again stable with the as:'mptotic period 
being finite as given by Eq.(197). The motion during this whole range of 
transition altitude is non-linear and a separate study should be undertaken 
to analyze all the Inherent phenomena. This would Include the frequency 
resonance and the eccentricity effect of the orbit. One simplifying aspect 
which significantly alleviates the difficulty in the analysis is that, at 
high altitude, the motion Is nearly Keplerlan. Hence, this leads to the 
pitch motion subject to gravity gradient and aerodynamic moments along an 
elliptic orbit. 

It should also be noticed that, if the orbit Is strictly circular, 
then from the linear equation (204), resonance effect is .lot detected. 

But, if the orbit is perturbed, it becomes slightly elliptic. In this case, 
by going back to Eq.(195) and expressing 6, 7, y and Ro/R as periodic func- 
tions In T, with period 2it, in the linearized form, we obtain a Mathieu's 
equation with a periodic forcing function having a period exactly equal to 
the pitch period. This will create resonance phenomenon which would 
ultimately build up larger amplitude for the perturbed angle-of-attack if 
the motion is uncontrolled. 
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Fig, 22. Limiting unstable altitude and resonance 
altitude in cruising flight. 

64 



Ballistic Entry 


Uni ike the case of cruising flight, in which the speed either remains 
constant or has a periodic variation with small amplitude, for ballistic or 
glide entry, effective entry is being made and as a consequence the speed 
decreases along the trajectory. In addition, the atmospheric density in- 
creases. This creates unsteady effects and, in the linearized form, the 
equation for the angle-of-attack oscillation is a linear second order 
differential equation with varying coefficients. 

We shall simplify the basic equation (195) to a form valid for both 
cases of ballistic and glide entry. First, we observe that the gravity 
gradient term in B and the following term in siny cosy are of the same order 
of magnitude. On the other hand, during effective entry, gravity gradient 
is small as compared to aerodynamic moment. Hence, with T = 0 and neglecting 
the two small terms mentioned, we are led to the equation 

“ C 6C 

a" +-^ a' + 6'C, + - 5(%) ^ cosy - = 0 (214) 

2v c y V 

Of the last three terms in this equation, because of the large factor 
the moment term is larger than the drag term. Therefore, we have the 
resulting equation 

C 

a" + ^ a' + 6'C, + 6C; - 5[(%) + 6C, C^] = 0 (215) 

2v *■ '- c y ^ 

Using the linearized expressions (200) and (202) for Cn» C(_ and we have 
the linear equation governing the variation of the angTe-of-attack 

a" + b(x)a' + c(x)a = - 6'C, + 6^C, C_ (216) 

L Q L o Uo 


where 


b(x) = 6[C, 


a 


(C + C ) 
m* m ' 

“ - 9 - 1 + 4 - 

2v 


T 


c(T) = - f (^) V - 

ot y C 


(217) 


a 


(X 


are two coefficients to be evaluated as functions of the independent variable 
X depending on the type of reentry trajectory. 


For the present case of ballistic entry, C, = 0. 
from Eq.(114) * 


On the other hand. 


Vl 

V 


2v 



(218) 
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Furthermore, from Eq.(n7) 


f- - f- - -IR. siny (219) 

where "ERo = -a » 900 for the Earth's atmosphere. Hence, from Eq.(217) we 
have 


b{x) = <5[C, - C. 

a ® 


— 

y 


c(t) • 6[C, a s1ny - {^) -r-2 - 6C. C, J 
'■(X c "y 


a 


(220) 


For ballistic entry at medium and large flight path angle, during the main 
portion of the reentry trajectory where the effect of deceleration is 
important, it has been shown in Chapter IV that the flight path angle remains 
nearly constant, that is siny = sin y©. Hence b(t) and c(t) are purely 
function of the dimensionless atmospheric density 5. This naturally leads 
to the change of independent variable from t to 5 through Eq.(219). Upon 
using this transformation we have the linear differential equation of the 
second order governing the angle-of-attack 

— + ^ + Cits)" = 0 (221) 


where 

bi = -^ + ki 
Cl = + ks 


with 


ki = 


1 


a siny, 


[C, - C, 


— >0 


k2 tC. 0 s1ny. - J > 0 

(a siny,)^ a c y 

k3 < 0 

(o sinyo)^ 


( 222 ) 


(223) 
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Concerning the signs of the constant coefficients k^. It Is trivial that ks 
1s negative. Since for entry a slnyo • -900 slnyo Is positive, the coefficient 
ki Is usually positive unless the damping in pitch Cmn, which Is always 
negative, Is negligibly small and the lift curve slope Clqj Is also small. 

For the coefficient k 2 , a negative Cma assures Its positiveness but an 
unstable, positive can make this coefficient negative. 

Since the linear differential equation (221) has varying coefficients, 
the positiveness of bi(6) and ci(6) Is not sufficient for the stability of 
the motion. To derive criteria for stability the recommended procedure Is 
to transform the equation Into a well known equation In mathematical physics. 

In this respect, Allen has considered a type of reentry missile such that 
upon some simplification the equation Is reduced to a Bessel's equation of 
order zero (ref. 18). On the other hand, Vinh and Laltone have transformed 
the equation Into a more general form of a confluent hypergeometric equation 
(ref. 19). In both these references the Independent variable used Is the 
altitude variable which necessitates the Introduction of exponential function 
Into the equation as an adequate representation for the atmospheric density. 
This inconven1en< e has been removed In the present formulation. Furthermore, 
we shall now use a general transformation of both Independent and dependent 
variables to globally discuss the stability of the angle-of-attack mode 
during ballistic entry for different aerodynamic characteristics. For this 
purpose, we consider the transformation 


a = e*^‘^y(n) 
6 = hn 


(224) 


where k and h are two constants to be selected such that the resulting 
equation has the desired form. Upon substituting Into Eq.(221), we obtain 


^ + [h(2k+k,) + 1] ^ + [h={k=+kk,*k,) 


. ]y . 0 

n ^ 


By selecting 



(225) 


(226) 


we have the equation 




= 0 


(227) 


If -Cfnoi Is not too small, the coefficient kj is much larger than ki and -kj. 
In this case, by neglecting the constant term in the coefficient of y in 
Eq.(227), and using the transformation 


n * 


(228) 


we are led to the equation 
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(229) 


^ + 1# + 2h(2kj 
dy2 U dy 


ki)y « 0 


By taking h such that 
2h(2k2 - ki) « 1 


(230) 


we have a Bessel's equation of order zero. Hence, for the type of ballistic 
entry bodies such that k 2 » ki, -ks, the approximate solution for the 
oscillation of the angle-of-attack Is 


a(6) - ■ (ki/2)]6 ) 

+ C2Y,(2 /tkj - (ki/2)]« ) ] 


(231) 


where ki and kz are given In Eq.(223) and Jo and Yo are the zero order 
Bessel's functions of the first and second kind, respectively. 

Since the Bessel's functions of order zero are oscillatory functions with 
decreasing amplitude and the density 5 Increases during ballistic entry, the 
condition for stability Is ki > 0, that Is 





(232) 


Going back to the more general equation (225), we now select the 
constants k and h such that 


h(2k + ki) » -1 
k^ + kik + ka = 0 


(233) 


The second equation In k has one positive and one negative root. We choose 
the negative root yielding 

k • - [1 ♦ /I - 4(k3/kn ] 

(234) 

^ - k. /T~. 4(kj/ki) 


The resulting linear equation for the angle-of-attack 1s a confluent hyper 
geometric equation 


n + (b-n) ^ - ay = 0 (235) 

dn" 


where b » 1 and 


a 


^ 2 /kl - 4k 3 


(236) 
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is a constant parameter, usually negative and large. The solution of (235) 
that will lead to a physically valid angle-of-attack variation as defined In 
Eq.(224) Is given by the Kummer's function 

y(n) ■ iFi(a, b, n) ■ 1 + ^ 

which In this case where b ■ 1, Is 

y(n) ■ 1 + an + n* + ... (238) 

( 2 !)' ( 3 !)^ 

Figure 23 shows the variation of the function y(n) for several values of the 
constant a. The function becomes more oscillatory for large negative value 
of a. From Eq.(236), and neglecting the small term k^, we have 

a - 1 - (k 2 /kx) (239) 


Hence, 1n general the ratio k 2 /ki Influences the oscillatory character of 
the angle-of-attack, while from Eq.(234) with k « -ki and, from the definition 
(224), the damping Is provided by the positiveness of the constant ki. We 
are led to the same condition (232) for stability. If we consider Eq.(238) 
as defining a family of curves depending on a parameter a, the best way to 
detect the oscillatory character of the function Is to select Integer values 
for -a » n, with n « 0, 1 , 2, ... Then we have 

yo - 1 

yi » 1 - n (240) 

y2 » 1 - 2n + 

« • • 

or in general 

Vn)- Jo ^ 


where 

n 1 n(n-l)(n-2)...(n-m-t-1) 

m m! 

\ / 


(242) 


Is the binomial coefficient. The solution Is a Laguerre polynomial, and as 
such, a higher order polynomial provides more zeros to the function. On the 
other hand, using h = 1/kx, we have n * ki6. The largest value of 6 
corresponds to sea level. Since the range of n Is large for large value of 
ki, the number of oscillations also Increases with ki. 

In summary for a general type of entry body the solution for the 
variation In the angle-of-attack Is given by 
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1 Fi Cl-(k2/ki ) , 1, ki6] 


(243) 


a “ a^e 


-ki6 


where ki and kz are given by Eq.(223) and the small effect of k) has been 
omitted. The fact that this solution tends to the Bessel's solution as given 
by Eq.(231) In the case where k2 » kj can be seen by using the asymptotic 
form for Kummer's function for large a 

iFi(a, b, n) - e^'^^ ^b-l^ /2bn-4an) (244) 


In the present case, b ■ 1, n"ki 6 , a ■ l-(k2/ki). Upon using this 
asymptotic approximation In Eq.(243), we have the solution (231). 


To verify numerically the analysis above, we Integrate the linear 
equation (221) using the following data. 


C» - 0.982 , 

Uo 


C, - 1.06 , • 0 

L m* 

a a 


C - -0.002 , C - -0.68 


m 


a 


m. 


(245) 


Yo • -20* , ky -1.129 . Ro/c -3.6498x10* 


This provides the basic numerical values 

ki - 0.0022 , k2 « 0.01026/ , kj - -1.099 x 10’^ (246) 

The Initial conditions are a® ■ 3®, a« * 0, at an Initial altitude such that 
6oC0p • 0.015. This value Is the initial deceleration In g's, due to the 
drag force at circular entry speed. This basic solution Is plotted In 
dashed line In Fig. 24. Several other solutions have been generated for 
different values of k2. While the damping per cycle remains essentially the 
same, the fact of Increasing k2, that Is the effect of the static derivative 
Cftu» produces more oscillations during the descent. This, In turn, shortens 
the period and leads to a stronger convergence In the angle-of-attack. For 
very large value of k2 the limiting curve tends to the Bessel's solution as 
first derived by Allen (ref. 18). 


From the definition (223) of k2, because of the large ratio Ro/c, this 
coefficient Is very sensitive to change In the static derivative 
Is to the location of the center-of- gravity. We have seen that k2 mainly 
governs the oscillatory character while the damping Is provided by the 
positiveness of the coefficient ki as expressed In condition (232). This In 
turn, depends on the damping In pitch which. In this case of the SSTO 
vehicle. In the absence of a horizontal Itablllze'r, must come entirely from 
the wing. 


To show that we still have convergence In the angle-of-attack while 
relaxing the condition on we consider two special cases with a positive 
value for First, we consider the case where k2 - k^. This corresponds 

to the value Cm^ - 3.64 x 10"‘*. Neglecting ka, we have a « 0 In Eq.(235). 
With b » 1, the exact solution of this equation Is 
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Fig. 24. Angle-of-attack oscillation durinq 
ballistic entry for various values of 
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( 247 ) 


y(n) ■ CiE1(n) Ca 

where E1(n) Is the exponential Integral 

E1(x) • Y + log X + J ^ 

n*l 

with Y • 0.57721... being the Euler's constant, 
n ■ 1 ( 1 $, the solution for the angle*of-attack Is 


(248) 


From Eq.(224), with k « -ki, 


a ■ [C.E)(k;J) ♦ Cz] (249) 

This, of course. Is a danced solution for ki > 0. 

The second case corresponds to ka » 0, which occurs for • 7.92 x lO'**. 
In this case a ■ 1 In Eq.(235) and we have the exact solution 


y(n) - e^ CCiEi(n) + Ca] (250) 

where the exponential Integral Is now 

Ei(x) • - Y - log x - I — (251) 

n*l 

The solution for the angle-of-attack Is 

a ■ CiEi(ki$) + Ca (252) 

Since the function Ei(x) Is decreasing, the pitch nwde is convergent but non 
oscillatory. 

Glide Entry 

The coefficients b and c and the forcing function In the basic linear 
differential equation (216) are now evaluated for the case of glide entry. 

As shown In Chapter IV, for equilibrium glide, we use the approximation 
Y * 0 and 


11 

V 





Hence, we deduce 



(253) 


(254) 


It appears that If the speed Is not too small, $' Is of the order of 5 and 
from Eq.(217), for the coefficient c(t) we can singly retain the moment term 
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because of the largt 'io Ro/c. This assumption Is cor.iistent with the 
assumption of equilib n glide since from the general equation (219) for 
the variation of the a« ity, when the speed is high, the flight path angle 
is negligioly small and 6* is small. Then, we write Eq.(216) 

a" + lc45a' + ks6a * -6'C, + 5^C, (255) 

Lo L.0 


where 


ki, 





) 

-] 


ks = 




(256) 


When is not negligibly small, the coefficient ks is large while the 
damping term k 4 remains small. By Eqs.(253) and (254) it is possible to 
change the independent variable to 5 as has been dene for the case of 
ballistic entry. But, in order to have a comparative analysis with the phugoid 
mode for g^ide entry obtained in Chapter IV, we shall use the kinetic energy 
V as the new independent variable. The transformation is done through the 
first of the equations (253) written as 

7' = . 11^ (257) 


where 


^ ° (258) 

is the lift-to-drag ratio used for the glide. Upon performing this change 
of independent variable in Eq.(255) and expressing 6 and 6' in terms of V, 
we obtain the linear differential equation of the second order governing the 
perturbation of the angle-of-attack during equilibrium glide 

7(1-7) ^ . [7 (l-V)] % + (259) 

dv“ ‘^'•0, dv ’'"0. 4v 

By the change of variable 

^ = 1 - V (260) 


with z increasing from zero for a starting glide at circular speed, we have 
the equation 


dl-;) — + [1 


- ( 1 - 


2C, 




a = 


E(2-0 

4TTCT 


(261) 


Identifying the homogeneous equation with the hypergeometric differential 
equation 
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C(l-d ^ + [c - (a+b+l);] ^ - aba . 0 

dr,2 


( 262 ) 


we see that 



ab = - ^ (263) 

c = 1 

A physically valid solution of the equation is given by Gauss series with 
c = 1 

F{a,b; t ... (264) 

(1!)^ (2!)" 

In general, both IC4 and ks are positive with ks » k^. Then a and b are 
nearly equal with opposite sign. Hence, w'th the approximation -b = a > 0, 
we have the solution for very large value of a 

F(a, -a; 1; = 1 - ? + -il- + ... (265) 

(1!)" (2!)^ 

= Jo(2a A) 


We recall that ^ varies in the range [0, 1]. The Bessel's solution is a 
damp ed soluti on and the number of oscillations increases proportionally with 
a = /Eks/4CDo. Hence, as in the case of the phugoid mode, higher lift-to- 
drag ratio provides more oscillations. Here, the main influence comes from 
a strong positive value of ks requiring a stable negative val’ie for 
For the case of glide entry, if the static stability derivative is relaxed 
thrusters reaction control must be supplemented to provide adequate stability. 

To display explicitly the damping and the oscillatory characteristics, 
we can perform a transformation of variables in Eq.(259), in a similar way 
as has been done in Eq.(184) f'r the phugoid mode. With the change of 
variable 

a = ■■■ -■ ■ ■ ^ , V = cos^U , ® S V* ^ ( 266 ) 

Vtanu/(cosy)*^‘*^^^o 

the linear equation (259) is transformed into 
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^ ! [{2 + s1n“u)‘ - (1 + 4 cos%)]}y = 

dy^ ^Do 4 sin^u cos^U 


E (1 •*• cos^y 
cos^y(cosy) 





( 267 ) 


If Is not negligibly small, ks is large. Then during the hypersonic 
portion of the glide where the speed ^s not small we can neglect the non- 
constant term In the coefficient nf y and obtain the solution to the 
homogeneous equation In the form 


1 

" 1 .. — 

V tany/(cosy 

To this equation, we add a particular solution of Eq.(261) to have the 
general solution. For small c we can consider the forcing function of this 
equation as a constant. Then it appears that the particular solution Is also 
a constant and Is very small 'for large kg. Hence, for all practical purposes, 
the angle-of-attack oscillation during equilibrium glide entry Is given by 
Eq.(268). The pitch mode is oscillatory and when v varies from 1 to 0, y 
varies from 0 to Tr/? and the number of oscillations Is approximately 


[Cl cos 




^ y) + C 2 sin (• 
Do 


Ek. 


y)] (268) 



(269) 


As for the damping. If k 4 > 0 the denominator of Eq.(268) Increases during 
the descent and the amplitude of oscillation of the angle-of-attack decreases. 
When k 4 < 0, this function increases If the only If 

1 + sin^y > 0 (270) 


Then, depending on the negative value of ki*, when y is large, that Is at low 
speed, the condition may not be verified and a can Increase after an Initial 
convergence. 

To verify numerically the analysis above, we choose a low angle-of-attack 
for glide with the following aerodynamic data 


= 0.301 , 

Lo 

= 0.452 , 

0 

II 

0 


= 0.186 , 


= 0.397 


(271) 


a 




=-0.012 , 


= -1.72 
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This provides the basic numerical values 

k4 * 1.6195 , ks = 3.879 X 103 ^ = 34.779 (272) 

The last number is the estimated number of oscillations if the glide is 
effected until sea level. Equation (261) is selected for the integration 
with the initial value a® = 3®, a® = 0 at ^ = 0.005 which corresponds to a 
slightly subcircular speed. The solution is plotted in Fig. 25 in dashed 
line. We notice the strong oscillatory character and danping of this pitch 
mode. They come mostly from the effect of Cma in the coefficient ks. As 
seen from condition (270), we still can have convergence of the solution 
with small and even negative value of k4. This means that if the location 
of the center-of-gravity is not far aft we can relax the effect of the 
damping in pitch C^q in the hypersonic glide mode. To test this assessment 
we have generated a^numerical solution, plotted in solid line in Fig. 25, 
with k4 = 0. According to the analytical solution (268) when ? h^s. increased 
fro m 0.0 05 to 0.4, which corresponds to a decrease of the speed (Y)^ = 

Vo/»^ gfto from 0.9975 to 0.7746 the amplitude of oscillation is reduced to 
0.2947 of the initial amplitude. This is in complete agreement with the 
numerical solution. 
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APPENDIX A 


CHARACTERISTICS F THE ATMOSPHERE 


In analytic formulation, there are two basic equations governing the 
density of the atmosphere as a function of the altitude. The first is the 
familiar equation of state for the atmosphere relating its pressure p, 
density p and temperature T 

p = pR*T (A.l) 

where R* is the gas constant (R* = 287.053 MKS for air). The second basic 
equation expresses that the rate of change of pressure must equal the increased 
weight of the atmosphere supported, as the altitude changes 

dp = -pg dH ^A.2) 


From the equation of state 


dp _ dp dT 
p ■ p ‘ T 

which, combined with Eqs.(A.l) and (A. 2), gives 


^ = .r -3- 

p >■ 


+ 1 HH 

R*T T dH 


The density gradient has been defined as 



or from Eq.(A.4) 


a = -R[ ] 


(A. 3) 


(A. 4) 


(A. 5) 


For the Earth, at altitude below 120 km, -a, which is a dimensionless 
quantity varies 'rom 750 to 1350 with a weighted average of 900. Fig. A.l 
plots the value of -a as function of the altitude. The discrete values are 
obtained from numerical differentiation of the tabulated U.S. Standard 
Atmosphere, 1962, while the smooth continuous curve is generated by using 
the approximate inverse polynomial representation of the atmospheric density 
as given in the U.S. Standard Atmosphere Supplements, 1966. It appears that 
the average value of -a = 900 as proposed by Chapman in the analysis of 
reentry trajectories can be used. 


At low altitude, the density gradient is used in the form 


1 dp a 1 r g 
p dH ■ R " ‘ T R* 



(A. 6) 


Using the standard atmosphere, in the troposphere, up to 11 km, dT/dH is 
equal to -6.5 x 10"^ degrees/m and a’'-,e 11 km, in the stratosphere, 
dT/dH » 0. Hence we consider the ,<ilues 
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H < n ,000 m , 
H > 11,000 m , 


1 dp _ 0.0276632 

p# T 

= -1 .57688 X 1 

p drl 


with T expressed In degrees Kelvin 


X 10 



Fig. A.l. The atmospheric density gradient as function of the altitude. 



APPENDIX B 


VEHICLE CHARACTERISTICS 

S » 6000 ft^ (557.42 m^) 

Z = 219.12 ft (66.79 m) 
c = 58.14 ft (17.72 m) 

W = 429,613 lbs 

I^ « 8.510,000 slug-ft^ 

ly - 50,970,000 slug-ft^ 

I^ » 48,600,000 s1ug-ft^ 

^xz 566,499 slug-ft^ 

Dimensionless Quantities 
B = 3(!j, - • -2.359623 

• 1.129 
R./c - 3.6498 X 105 

Aerodynamic Characteristics 

As provided in reference 8, with typical values given in 
Eqs.(245) and (271). 
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